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ABSTRACT 


A new two-dimensional data-adaptive algorithm utilizing the iterative Toeplitz 
approximation method is presented. This algorithm provides a robust and efficient 
means for accurate estimation of 2~D autoregressive parameters representing spatially 
variant data arrays. Its convergence performance is comparable to that of the 2—-D 
Recursive Least Squares (RLS) algorithm but is computationally more efficient. The 
savings in computation is realized by reducing the size of the matrix to be inverted 
when solving the AR model normal equations. The ability of the algorithm to accu- 
rately estimate the model parameters using very small data sets and various window- 
ing schemes are evaluated. Spectral estimates are compared for quarter-plane (QP), 
nonsymmetric half-plane (NSHP) and combined-quadrant (CQ) regions of support. 
Additionally, the algorithm is tested in noise cancellation and line enhancement ap- 
plications using image arrays. This algorithm may be implemented for data-adaptive 


image processing or coding and for applications requiring 2-D spectral estimation. 
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I. INTRODUCTION 


A. PARAMETER ESTIMATION OF DATA-ADAPTIVE ARRAYS 
Two-dimensional (2-D) autoregressive (AR) modeling provides a means for the 
spectral estimation of signals received by an array of spatially distributed sensors. 
Additionally, the 2-D AR model parameters may be used to represent the original 
signal in a compact manner resulting in compression of the data to be processed or 
transmitted. This is particularly useful in many image processing applications. The 
traditional least squares or Wiener filtering based 2-D AR parameter estimation [Ref. 
1] provide accurate spectral estimates but require the the inversion of a relatively large 
autocorrelation matrix. This is undesirable since matrix inversion is computationally 
expensive. Iterative methods, using a Toeplitz approximation algorithm, have been 
suggested (Ref. 2, 3] as an alternative to the direct Wiener filter solution. The itera- 
tive Toeplitz approximation method has proven to be an effective means to estimate 
the parameters of fixed data arrays without having to invert the correlation matrix. 
This method is particularly useful when data arrays of limited size must be processed. 
The results obtained using iterative methods to obtain the AR parameters of 
fixed data arrays suggest that these methods may be used to some advantage when 
estimating parameters adaptively. Adaptive parameter estimation is necessary in 
many digital signal processing applications where the data is non-stationary. In these 
applications, real-time or near real-time processing is often desirable. This requires 
that any algorithm used must be computationally efficient and be capable of providing 
accurate estimates obtained from a minimum of data. This thesis addresses the 


implementation of the iterative Toeplitz approximation method for 2-D parameter 


estimation of non-stationary data. It is shown that this algorithm can be efficiently 


implemented to obtain accurate estimates from very small data arrays. 


B. EXPLANATION OF NOTATION 

All rectangular matrices are denoted by a boldface, captial letter, e.g., R. Col- 
umn vectors are designated by a boldface lowercase letter, for example: x. Occasion- 
ally, a vector created temporarily for the development of a mathematical expression 
will be represented by a lowercase, boldface, greek letter. An example of this would 
be the vector a. A special operator matrix will be represented by a calligraphic cap- 
ital letter, such as F. Scalar values are generally represented non-boldface lowercase 
arabic or greek letters, e.g., a or A. Captital letters such as A or P are used to rep- 
resent a variety of values or expressions. These include anything from a polynomial 
to the dimension of a matrix. The symbol T' is reserved to indicate the transpose of 


a vector or matrix. 


C. THESIS OUTLINE 

The following describes the organization of the remainder this thesis. Chapter 
II serves as an introduction to the problem of 2-D AR modeling and parameter esti- 
mation. Particular attention is given to the formulation of the normal equations for 
quarter-plane (QP) and nonsymmetric half-plane (NSHP) regions of support. This 
chapter concludes with a description of the least squares algorithm and direct in- 
version method of solving the normal equations. This serves as the foundation for 
Chapter III which extends the concept of solving normal equations using the itera- 
tive Toeplitz approximation method to derive parameter estimates referred to as the 
fixed data method. Chapter IV proceeds with the development of the data-adaptive 
algorithm using Toeplitz approximation. Two-dimensional AR spectral estimation 


is discussed in Chapter V, along with a comparison of results using both fixed and 


adaptive iterative methods for spectral estimation. Additionally, a discussion of the 
use of the data-adaptive algorithm for noise cancellation and line enhancement of 
image arrays is provided in Chapter VI, with experimental results. Conclusions and 


recommendations for future work are found in Chapter VII. 


Il. 2-D AUTOREGRESSIVE MODELING 


A. OVERVIEW 
This procedure assumes a stationary (homogeneous) random process z[nj, 79] 
that is the response of an AR model excited by a white noise input w[n,,n2| having 


a variance a2. The AR model as shown in Figure 2.1 may be described as an all pole 


e 


w(n,, 72) Homes) r(n1,N2) 


Figure 2.1: AR model excited by white noise. 


filter with the transfer function 


1 
1¢+> aan ko €A) a(ky, kz )e—3 (wr ks tw2ko) 


where A is the region of support over which the parameters a(k,,k2) are non-zero. 


H(w,,w2) = (25) 


The difference equation for the system that generates r[n,, 2] can be written as 
z[ni,n2] = — >_> a(2,7)2[ny — 2, ng — J] + w[mi, no]. (222) 
(i,j) €A 
A linear set of equations for the filter coefficients a(z, 7) may be formed by multiplying 
both sides of (2.2) by z[n; —1,,n2 — 12] and computing the the statistical expectation 
of the resulting expression [Ref. 4]. This leads to the following expression called the 
normal equation 


Re{h, l,] i ys S5a(7, 7) Re[h = 2, ly ae J ’ (2.3) 
(t,7) EA 


for 1;,l2 > 0. The coefficients a(z,7) can be derived from this normal equation. It 
should be noted that the structure of the normal equation depends on the ore of 
support A. A region of support may have any shape, but thé most commonly used 
are the quarter-plane and nonsymmetric half-plane. 
1. Quarter-Plane Support 
The most straightforward region of support is the quarter-plane. A region 


A is considered to have quarter-plane (QP) support when a(z,7) has non-zero values 


in one quadrant only as shown in Figure 2.2. For this case the normal equation has 





Figure 2.2: The quarter-plane region of support. 


the form 
Py-1 P2-1 


Rell) = S> So a(t, 7) Refl; — i,l2 — J], (i, 7) ¥ (0,0) (2.4) 


1=0 7=0 


where /, = 1,2,...,—land l, = 1,2,..., P2—1 with P; and P, being the dimensions 


of A. If it is assumed that a(0,0) = 1, we may express (2.4) in block matriz form as 


Ro R_; R_» eee R_p,41 ao € (9) 

R, Ro R_; cbt R_p, +2 ay 0 

R» Ry Ro --- R_pys a2 |=| 0 |. (2.5) 
Rp,-1 Rp,-2 Rp, =3 "* Ro AP, -1 0 


Each block Ray of this matrix is given by 


R(M,0) R(M,-1) +: R(M,—P) +1) 
Ruy = RP, = RUM, 1) RM, 0) 7 ae! ai +2) on 
R(M,P;-1) R(M,P,—2) --- — R(M,0) | 


while the coefficients are ayy = [a(M,0) a(M,1) a(M,2)...a(M, P. — 1)]? and the , 
error variance vector is €(% = [o?0...0]7. The blocks Ry along the diagonals of the 
autocorrelation matrix R are equal and the diagonal elements of Ry are equal, thus | 
R is Toeplitz block- Toeplitz. 
2. Nonsymmetric Half-Plane Support 
A region A with non-zero a(z,7) as shown in Figure 2.3 is considered to 


have nonsymmetric half-plane (NSHP) support. The normal equation for this case 





Figure 2.3: The nonsymmetric half-plane region of support. 


must be modified to account for non-zero a(0,7) in the fourth quadrant. This may 


be written as 


b2+FP2-1 P,-1 P2-1 
Relh,b}= © a(0,3) Rell, ls —3] + DE De alg) Re[4 eats — 9), ee 
j=1 20 3=0 


minete(@ gi Ue — Ol 2 acne | and 


oe 0,1,2,...,.fb2.4+P-1, for /; = 0 
eee eo... bo tort, 0 


The dimensions of A are given by P, and P2 and Lz is a negative number as defined 
in Figure 2.3. As in the quarter-plane case, if we let a(0,0) = 1 we may write (2.7) 


in block matrix form as 


Ro Ri R22 -: Repay aR &(0) 

R, Ro R_; ses Roop 4e a) 0 

R2 Ri Ro oe ops a = 0 (2.8) 
Rp,-1 Rp, —2 Rp, —-3 =6°"" Ro aP,—1 0 


The blocks Ro, Ry and Ry each have different forms. The diagonal block 
R, has the form 


R(0, 0) R(0, -1) --- R(O,—-L2 —- P2 +1) 
s R(0, 1 R(0, 0 “> RO, -L2 — P2 +2 
R, = . ) ) ; ( 2 2 + 2) (29) 
R(0,l2+P2-1) R(O,L2+P2-2) --- R(0, 0) 
the blocks Ray = R7,, in the first row and column may be written as 
R(M, —L2) R(M, —L2 -1) --- R(M,—L2 — P. +1) 
R(M,-L2+4+1) R(M, —L2) --» R(M,—L2 — Po +2) 
. ; . , (2.10) 
R(M,-L.+ P2-1) R(M,-I.+P,-2) --- R(M,0) 


and the remaining blocks Ry are given by (2.6). 


The model parameter vectors are given by 
Ao = [a(0,0) a(0, 1) a(0,2)...a(0, LZ. + Pp —1)]? 


and 


am = [a(M, L2) a(M, Lo +1) a(M, L. + 2)...a(M,L, + P, —1)]7 


for M = 1,2,...,P,—1. The error variance vector €) = [o?0...0]7. Except for the 
upper diagonal block (2.9) and top and left most blocks (2.10), the NSHP autocorre- 
lation matrix is block-toeplitz with each block being toeplitz as well. Quarter-Plane 


support can be considered a special case of NSHP support with Lz = 0. 


B. SOLVING THE NORMAL EQUATIONS BY DIRECT INVERSION 
The normal equations also arise in the 2-D linear prediction problem. When 
solving the normal equations, the objective is to find the parameters which mini- 
mize the prediction error [Ref. 5]. It is frequently more convenient to describe the 
formulation of the normal equations in terms of linear prediction. 
In the 2—D least squares problem, the error between the random process z[nj, n2] 
and its estimate Z[n, n| is given by 
e[ny, no] = z[n1, n2] — t[n1, no] (2.11) 
or in expanded form as 
e[ni, n2] = z[n1, 72] aie ~ » a(2,7)z[ny a 2, N2 oa a (2,7) # (0, 0). (2.12) 
(1,7) €A 
The objective of the least squares method is to minimize the squared errors from a 
particular set of these terms [Ref. 4]. If we let a(0,0) = 1, we can express (2.12) as 
e[ny, n2| = > » a(z,7)z[ny = 2, nga—- J (213) 
(i,3) €A 
which is compactly represented in vector notation as e = Xa. The error e[nj, 72] 
is computed for each position of the filter, then normal equations are formed by 
multiplying both sides of (2.13) by X7 and applying the orthogonality principle [Ref. 
5]. This results in 


Ra=€ (2.14) 


where R = X7X and the error variance vector € = X7e = [€0)7,0,...,0]’ with 
€) = [a?,0,...,0]7. The matrix X may be defined as 


|| | 
>. e— Xo > a Xo Xi °° XP,-1 (2.15) 
| 


and the model parameters may be given bya=[1 a’ ie The parameter estimates 
are obtained from (2.14) by muliplying both sides of the equation by the inverse of 


R which may be written in expanded form as 
A= (XM ee Xe (2.16) 


This is referred to as the direct inversion method. 

The preceding discussion assumes that R is Toeplitz block-Toeplitz. This true 
for the case where the autocorrelation method (Ref. 5] is used to form the correlation 
matrix or when the theoretical correlation is known. In practicality, however, param- 
eter estimates must be obtained from relatively small sets of finite data. In this case, 
it is often desirable to compute R using the covariance method. This results in a 
block autocorrelation matrix R with non-Toeplitz blocks Ray [Ref. 2]. The rows of 
the matrix X are the elements of the random process covered by the filter mask for 
each point being estimated as the mask is moved over the data. In the covariance 
method the filter mask is not moved past the boundaries of the region of support. 
This means that when QP support is used, X will be of dimension P, P) x P,P. where 
P, and FP, are the dimensions of the data array accessible to the filter mask. For 
example, using a 3 x 3 (nine element) filter mask to form the normal equations of an 
8 x § data array would result in an X matrix with dimension 36 x 9. It should be 
noted that P; = P, = 6, since two rows and two columns of data cannot be estimated 


by the filter. 


The preceding discussion provides the basis for the subsequent work in this 


9 


thesis. The next chapter will addresses this same problem using the iterative Toeplitz 


approximation method. 
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Ill. ITERATIVE TOEPLITZ 
APPROXIMATION ALGORITHM 


The major drawback of the least squares method described in the previous 
chapter comes from the requirement to invert a relatively large autocorrelation ma- 
trix. This is undesireable since matrix inversion is computationally expensive. In 
this chapter we present an iterative method for solving the normal equations which 
does not involve the direct inversion of the correlation matrix to obtain the model 
parameters. This method is then extended to the case where the covariance method 
is used to form the correlation matrix and the normal equations are solved iteratively 


using Toeplitz approximation. 


A. THE ITERATIVE SOLUTION OF THE NORMAL EQUATIONS 

An alternative to direct inversion is to take advantage of the near Toeplitz- 
block-Toeplitz structure of the autocorrelation matrix and successively partition the 
normal equation [Ref. 2]. This partitioning permits the estimation of parameters 
while requiring the inversion of a matrix that is significantly smaller than the original 
autocorrelation matrix. The following discussion develops the iterative solution for 
QP and NSHP regions of support. 

1. Quarter-Plane Support 

In order to iteratively solve the QP normal equations we must begin by 


dividing both sides of (2.5) by o?. This results in a modified normal equation given 


a 


by 


Ro Re, Reo eter a ag € (0) 

R, Ro R=) ---eRSpe> ay 0 

R» R, Ro «+: R-p, 4s ay |=| 0 |. (3.1) 
Rp) Rpew Reason | ets 0 


where €) = [100,...,0]? and a’ = aj for2 = 1,2,..., PR; —1. We now partition 


the normal equation (3.1) as follows 


nf Rel lees! = Lo" 
Pe; Airy Oe (3.2) 
where 
Ro Rope Re Reps 
R; Ro R.,  --- R pis 
G,=| R R, Ro +--+ Reps] , (3.3) 
Rp, -2 Rp, -3 Rp, -4 ae Ro 
hp =(Rp-1 Rp». Res --- Ril, (3.4) 
a, =(a? ae! |. (3.5) 
and 
=(€OT oF ... OF)’. (3.6) 


We may now form a set of coupled iterative equations by defining an op- 


erator F, which is given as 


Gr’ OO 
Fi=| 0 a ’ (3.7) 
and using this operator to premultiply both sides of (3.2) to yield 
a1=Gr'fy, —hn ap, —il (3.8) 
and 
ap, 1 = —Ro'hjar. (3.9) 


2 


Equations (3.8) and (3.9) suggest an iterative solution and can be written as 


oe Gala heen), (3.10) 


and 


an”), = —Rpthlal | (3.11) 


which provide a means to solve the normal equation (3.1) iteratively. 

It can be noted that the submatrix G, is nearly the same dimension as 
the original correlation matrix. This means that the inversion of G, provides little 
advantage, computationally, over the direct inversion method. It is clear, however, 
that further partitioning of G, will decrease the required complexity. Partitioning of 


G, yields 


nf Roles’) =['o' 
aes eels | (3.12) 
where 
Ro R_; R_2 --- R_p4s 
R, Ro R_, +--+ Rop4s 
G, = R» R, Ro lie R_p, +5 , (3.13) 
Rp,-3 Rp-4 Rp,-s -:: Ro 
= | Rp, -2 Rp, -3 Rp, -4 eee R; | 9 (3.14) 
ana al? Ay) Pes Sarees ae (el5)) 
and 
ea nO ee Oe (3.16) 


As before, both sides of (3.12) are premultiplied by an operator defined by 


G;' a 


F, = |%3 R>? 


(3.17) 
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which results in 


a2 = Gy'[y 2 — he ap,_2], (3.18) 
and 
Ap, 2 => —Rj"hi a2 : (3.19) 


Equations (3.18) and (3.19) may be solved iteratively in the same manner as (3.8) 


and (3.9). This solution is expressed as 
a$? = Gp"fy, —heap®s”], (3.20) 


and 
ap = ee (3.21) 


This partitioning process may be continued until Gp,_1 = Ro. This will 


result in a partitioned normal equation given by 


Gp,-1 a oe “oo 
! i : = 1 gee 
| hp, -1 Ro ap_» 0 ( ) 


where hf _, = Ri, ap-1 =agandy p_, = E. 
Premultiplication of both sides of (3.22) by an operator Fp,_, will result 


in the iterative solution 


aj” = Ro —[(Rja,"" +... + Rees 2 |, (3.23) 
and 
at) = —R5{Riag 9) + Raay) + ...4+ R_pscans”)- (3.24) 


For this case the operator matrix is given by 


Paso Ra" el , (3.25) 
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At this point in our development, we can combine the preceding iterative 
solutions from the successive levels of partitioning into a set of compact iterative 


summations 


all(*) — a) 2 Re x R_,a’ Wk- te (3.26) 


af = RSS Reval 45), (3.27) 


+=0 


where af) = Ro'E), a) = —Ro'R;ap, and j =1,2,...,P,) — 1. The index of 
iteration is k and a’ are the P, x 1 parameter vectors. The solution of (3.26)-(3.27) 
requires O(P?P?K) multiplications, where K represents the number of iterations 
to converge to the true parameter values. An additional P,P? multiplications are 
required to compute the initial parameter estimates ae fOr ey =", 2h ieee 
This results in total multplications of O( P?P?K + P,P?). Depending on the number 
of iterations required to converge to the true parameters, this total compares with the 
algorithms of Wax and Kalaith [Ref. 6] and Akaike [Ref. 7] which require O(P? P?) 
operations. When large arrays are used this represents a noticeable savings over direct 
inversion which requires O( P? P;) multiplications. 
2. Nonsymmetric Half-Plane Support 
The derivation of the iterative NSHP solution follows a procedure of suc- 
cesive partitioning similar to that of QP support. However some differences exist due 
to the asymmetry of the autocorrelation matrix consisting of three different types 
of block matrices. This asymmetry results in a somewhat modified set of recursive 


summations given by 


al) — aio) =H — R_ iallk- -1) (3.28) 


t=] 


P;-1 
ai(*) = =Re Ras! = Re s~ R,_,a(*-) | (i = 5), (3.29) 


t=1 
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where ap) = Rp'EO, ae = —R5'R, a0, and j7 = 1,2,...,4% —1. The mdex 
of iteration is k and a’! are the P2 x 1 parameter vectors. As with QP support, The 
solution of (3.28)-(3.29) requires O(P?P}K) multiplications, where K represents the 


number of iterations to converge to the true parameter values. 


B. THE ITERATIVE SOLUTION USING TOEPLITZ APPROXIMA- 

TION 

The discussion of the previous section assumes an autocorrelation matrix with 
Toeplitz-block-Toeplitz structure. In most cases the autocorrelation matrix must be 
formed from a limited, and possibly very small, amount of sampled data. In this 
situation, it is best to use the covariance method to estimate the correlation matrix 
to reduce the bias in the estimate. However, since this method does not have Toeplitz- 
block-Toeplitz structure, it is necessary to modify the iterative algorithm described 
above. 

1. Forward Iteration 


For first quadrant QP support, the covariance estimate R has the form 


Roo Ro, Ro2 +: Rop-i 

Rio Ri 2 ee 

R20 R21 Roo 3 --: Ro p-1 (3.30) 
Rp-io Reais Rp-i2 --- Re-i,p-i 


A proven procedure [Ref. 8] that can be used to form the Toeplitz approx- 
imation is to first average the diagonal blocks R;;, 2 = 7 and then form a Toeplitz 
T matrix from this averaged block R,,, by averaging its diagonal elements. The 


diagonal elements of T can be obtained as follows 


: Pig 
aa > Ravg(? + 3,7) (3.31) 


j=0 





where 2 = 0,1,...,)—1. 
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We may now proceed to succesively partition the normal equations in the 


manner outlined in the previous section. The first partitioning leads to the iterative 


equations 

a= Gay hua |e (3.32) 
and 

ape = —ReLpahper (3.33) 


The diagonal block Rp, -1,p,-1 may be written as 
Lee ee A ne (3.34) 


where Ap, -1 p,-1 is the difference between Rp, -1,p,-; and the Toeplitz approximation 


fl be 
Substituting (3.34) into (3.33) modifies the recursion to give 
a}? = Gi'fy , —hi apt”), (3.35) 
and 
ap), = —-T Ap-17-1an'.) — Thrall). (3.36) 


The final expressions resulting from the successive partitioning are then 


k = k= 
pe Ge ype hae ase (3.37) 


and 


al(*) = —R,ahp-iae-)) ‘ (3.38) 


where ap,-; = a9,Gp,-1 = Roo, and hp,-1 = Roi. The diagonals R;,; may be 


rewritten as 


R;; = T+ A;,;, (3.39) 
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Now, substituting (3.39) into (3.38) results in 


= k= 
av, = Gry leper ea anele (3.40) 
and 
pie mens UVa ei) Td Srey Ae (3.41) 


Finally, combining the succesive iterative equations results in our solution to the 


normal equations using the forward iterative algorithm for the QP case, which takes 


the form 
P;-1 
an) = a) a T- Annan _7T7 3 Ry,a,0-)) (3.42) 
t=1 
i(k) T-1A (k—1) T-1 ket, R "(k-1) 3 
AS 5,59; = > Ria; (3.43) 
s=0 
tj 


where ag = T-1¢(), ey = TR;oag° and j = 1,2,...,%—1. The forward 
iteration will provide parameters for the first quadrant spectral estimate. 
2. Backward Iteration 

A filter mask may have quarter-plane support in any two quadrants. A sec- 
ond quadrant estimate exists for the quadrant adjacent to either side of the quadrant 
designated to be the first quadrant. The Hermitian symmetry and Toeplitz-block- 
Toeplitz nature of the autocorrelation matrix causes the estimates produced from 
diagonally adjacent quadrants to be identical [Ref. 9]. 

The second quadrant AR parameters 6(2, 7) may be estimated from the nor- 
mal equations using the backward iterative Toeplitz approximation. When 0(0,0) = 1 


the second quadrant normal equation is given by 


Roo Roi Ro sie Ro,P,-1 bp,-1 0 

Rio Ri Ri2 mes Ri P,-1 bp,-2 0 

R20 R21 Roa2  -*: | Rape bp,-3|=] O |. (3.44) 
Rp,-10 Rp-ia Rei -:: Rep-1,p,-1 bo E() 
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where R;; = R7,, €© = {100,...,0]? and 


jst? 


b? = (1/o7)[b(k, 0), b(k, 1), b(k, 2),...,(k, Py —1)) . 


It may be noted that the second quadrant autocorrelation matrix is identical to that 
of the first quadrant. However, the parameter estimates (7,7) are not generally the 
same. The estimation of the backward parameters is computed iteratively in the same 
manner as the forward parameters, but, the backward method differs in the way that 
the normal equation is partitioned. In this case the partitioning begins at the upper 
left diagonal block and continues until the G matrix consists of only the lower right 
diagonal block. The Toeplitz approximation T is identical to that of first quadrant 
support. 

The combined backward iterative equations, which are similar to those for 


the forward iteration, can be written as 


P,-1 
by” = by? — T7"Ap-1a-iby? — TT? > Ra-ap-1-bS 2), (3.45) 
7=1 
k ea — (k-1 
bY) = —T Ap 1-5,7R,-1-jbf) — T7? SY Ra-i-j,p,-a-sb G7), (3.46) 
s=0 
if) 


where by? = T71€, b® = T-Rp_yjb©) and j = 1,2,...,P;—1. The first and 
second quadrant parameters are used to find the combined-quadrant spectral estimate, 
which is discussed in Chapter V. 
3. Iterative Method for NSHP Support 
The iterative Toeplitz approximation method may be used to estimate the 
parameters of arrays with nonsymmetric half-plane support. As shown previously the 


NSHP autocorrelation matrix has an asymmetric structure which can be expressed 


Ibs) 


Ro,o Row Roo --- Ro p-1 

Rio Ri ioe eet 

R20 R21 Ro2 -+: Rope (3.47) 
Rp,-10 Rp,-1,1 Rp, -1,2 wot Rp, -1,P,-1 


The dimensions of the top most and left most blocks are reduced by L2 rows and 
columns, respectively. As in the previous cases the iterative process begins with 
forming the Toeplitz approximation of Ray. In the NSHP case, however, the upper 
left main diagonal block is not included in the average becuase it does not have 
the same dimension as the other main diagonal blocks. Thus, it is necessary to use 
the elements of the upper left main diagonal block to create a separate Toeplitz 
approximation designated as As 

The normal equation is successively partitioned as described in the previous 
sections. At each stage of the partitioning, the blocks along the main diagonal are 
expressed as the sum of the Toeplitz approximation and a difference matrix. The 


iterative summations are given by 


é - yi P,-1 
at*) = a0) _f Mogaea open > Roane (3.48) 
i=! 
is a P,-1 
al) TA. a9) Sas eee | (SBR ame (3.49) 
+=0 
tj 


where an — T71E(0), ARY = TR; a0” and 3 = 1,2... 


C. SUMMARY 

It is clear that iterative methods have a computational advantage over the di- 
rect inversion method. This has been shown for the Toeplitz-block-Toeplitz case and 
extended to the case where Toeplitz approximation is used to compensate for the 


non- Toeplitz nature of autocorrelation matrices formed using the covariance method. 
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Generally, the matrix that must be inverted using the iterative method has the di- 
mension of a single block in the block matrix. This corresponds to the filter order or 
mask size. The mask size will vary with the order necessary for accurate AR mod- 
eling of the sampled data, but filter masks that range from 3 x 3 to 6 x 6 will be 
sufficient for most applications. Although other methods may exhibit faster conver- 
gence to the true parameters, they do so with greater computational complexity. The 
iterative methods provide an approximation much sooner. Additionally, the iterative 
methods are guaranteed to converge when the block matrix is symmetric and positive 
definite [Ref. 2]. This property is important when extending the iterative Toeplitz 
approximation algorithm to the data-adaptive case which is introduced in the next 


chapter. 
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IV. DATA-ADAPTIVE ITERATIVE TOEPLITZ 
APPROXIMATION ALGORITHM 


A. OVERVIEW 

The necessity for adaptive filtering techniques arises when it is desired to process 
signals that result from an environment of unknown statistics [Ref. 10]. There exists a 
broad class of problems that fall into this category, these include such diverse fields as 
sonar, radar, image processing, seismology and communications. In general, adaptive 
filters provide a significant improvement in performance over fixed coefficient filters 
designed to operate on data with unknown statistics. Additionally, the use of adaptive 
filters makes possible new signal processing capabilties that would not be available 
otherwise. 

One distinct advantage associated with adaptive filters pertains to their ability 
process data on-line. This is desirable for many applications such as autoregressive 
spectrum analysis, detection of a signal in noise, adaptive noise cancellation and line 
enhancement. 

Two-dimensional adaptive filters are used to process data obtained from spa- 
tially variant data arrays. The most successful algorithms currently implemented for 
this purpose include the 2-D least mean square (LMS) and recursive least squares 
(RLS) algorithms. Both of these algorithms have been derived from the 2~D Wiener 
filter and method of least squares [Ref. 11]. With this in mind, it follows that the 
iterative method for fixed data may be successfully extended to operate on spatially 
variant arrays. The remainder of this chapter describes the development of the data- 
adaptive iterative Toeplitz approximation algorithms and some considerations that 


apply when using this algorithm. 


oe 


B. THE DATA-ADAPTIVE ALGORITHM 

The data-adaptive Toeplitz approximation algorithm is based on the same suc- 
cesive partitioning described for the fixed data case. The algorithm becomes adaptive 
in nature when the autocorrelation matrix R, is updated for each shift n of the filter 
mask. This yields a continuously updated set of AR parameter estimates. Compu- 
tation of AR parameter estimates begins with the first iteration [Ref. 12]. As in the 
fixed data case, only the approximated Toeplitz matrix T is inverted. Computation 
time is shortened further by the fact that R, is computed only from data covered 
by the filter mask, which amounts to taking the outer product of a P,P: x 1 vector 
where P,P: is equal to the number of filter coefficients. The data-adaptive form of 


the combined iterative equations can be written as 


P,-1 
al) = al) _ d him [Sagan ee as aloe S~ Revgiar a! (4.1) 
1=1 
P,-1 
n ~- n—1 - (n—1 
a Aaa oe Ran (4.2) 
1=0 
if; 


where aii) = T-1¢(°) | ere = Moartarnag andj = 1,2,...,P,—1l. These can be 
compared to (3.42) and (3.43) for the fixed data case. The initial parameter estimates 
A are determined for the initial point being estimated and all subsequent parameters 
a’ are a function of the parameters computed at the previous mask positions. The 


backward equations, used for the second quadrant estimates, can be written as 


P,-1 
bi” = bj? — Tz*A,p-1,7-1by” ? — Tz’ > Ran-taa-ibeiy;, (4.3) 
t=1 


P,-1 

be = —T 71 Ay,P,-1-5,,-1-3b0" ~~ sli 2 R,.P,-1-5,P; oe Dae (4.4) 
1=0 
ix) 


where b§? = TZ1E), bo!) = To7Rop,_1,;b?) and j = 1,2,...,P, —1. 
It should be noted that the only difference between the data-adaptive equations 


and the fixed data equations is the the index of iteration which is n instead of k 
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and the manner in which the autocorrelation matrix is formulated. In the data 
adaptive case, the iteration proceeds with each shift of the filter mask or update of 


the autocorrelation matrix which has the recursive form 
R, = Ra-1 + XnX;, - (4.5) 


The P,P. x 1 vector x, is obtained by arranging the 2-D data covered by the filter 
mask at each shift n = (n;,n2) in a vector form. The diagonal blocks of R, are used 
to form T,, in the same manner as the fixed data case. 

As with the RLS algorithm and other recursive implementations of the method 
of least squares, it is necessary to introduce a weighting or forgetting factor when 
computing R, [Ref. 10]. The forgetting factor is a scalar value that may be designated 
as 6(n) where n is the iteration number of the point being estimated. The weighting 


factor 8(n) has the property that 
0< A(n) <1, n=1,2,...,p, (4.6) 


where p is the total number of points being estimated. The purpose of {(n) is to 
ensure that data in the aietant past is ‘forgotten’ which will make it possible for 
the filter to adapt to statistical variations of the observed data when operating in a 
non-homogeneous environment. A commonly used form of the weighting factor is the 


exponential weighting factor which is defined as 
1d lp a aims (Aro Fi: (4.7) 


where 4 is positive constant with a value less than 1. The quantity 1/(1 — A) can be 
considered as an approximate measure of the memory of the algorithm. When 4 = 1 
we have the ordinary method of least squares which corresponds to infinite memory. 


By this reasoning, we may introduce a method of exponentially weighted least squares, 
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where {(n) is used as a weighting factor in the expression of the performance measure 


(Ref. 10] 


p 
E(n) = DU AP Ie(n)/* , (4.8) 
n=1 
where e(n) is the error defined by (2.12) at mask position (ni, 72). 
To see how the forgetting factor is introduced in the recursion we must refer to 


the original definition of the normal equation 
eae , (4.9) 


where a is the vector of optimum parameters which results in the minimum value of 
the performance measure (n). The P,P, x P,P, autocorrelation matrix in this case 
is defined by 
R, = 3 p\ ai ap ae (4.10) 
n=1 
We may now isolate the term of (4.10) that corresponds to n = p from the rest of the 


summation on the right side of the expression to obtain 


p-1 
R, = b dase ab cai cae + XnxX2 : (4.11) 


n=1 





The term inside the brackets on the right side of (4.11) is actually the sum of the pre- 
viously computed autocorrelation matrices R(n — 1), therefore we have the recursive 
relationship 


R,, = ARn-1 + Xnx2 . (4.12) 


The affect of 6(n) is apparent when (4.12) is expanded to give 


R, = xix) 

R, = Axixt + Xx? 

R3 = ?x xi + Axoxt + x3x? 

BN ey Oe xe x (4.13) 
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which clearly shows the decreasing contribution made by past data vectors when 
Areal: 

The choice of a value for is dependent on the statistical nature of the data 
being processed. Generally, if the data is random or uncorrelated it is undesirable to 
use the past data to compute the current estimate, therefore a small value for \ would 
be the appropriate choice. Conversely, a value of \ close to 1 would be desirable for 
highly correlated data. 

Since the statistical nature of the data is often not known a prior: when using 
data-adaptive filters, it is difficult to choose an optimum value for the forgetting 
factor. A further modification of the recursive equation (4.12) has been implemented 
for the data-adaptive iterative algorithm to address this problem. This modification 
involves the application of a weighting factor to the current update of R,. This 
weighting factor depends on the index of iteration n such that its value approaches 
1 as n —> oo. This biasing scheme causes the influence of the present data on the 
parameter estimates to increase with the number of iterations. The modified recursion 


is given by 





—] 
R, = R,-1 + (7 ; ) xnxl (4.14) 
The weighting factor on the second term in (4.14) is added to provide a more gradual 


shift of influence to recent updates than would be possible with only the » term 


present. This factor is not as important when a highly correlated signal is processed. 


C. SUMMARY 

The adaptive Toeplitz approximation algorithm can be used to estimate param- 
eters for the same regions of support as in the fixed data case. Additionally, various 
schemes for pre-windowing or post-windowing the data may be employed to improve 


the estimates. Other factors affecting the performance of the algorithm include the 
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directional scheme for moving the filter over the data array and the criteria used to 
determine convergence of the parameter estimates. This algorithm has been used for 
spectral estimation, image noise cancellation, and line enhancement. Experimental 


results involving these applications are presented in the next two chapters. 
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V. 2-D AR SPECTRAL ESTIMATION 


A. OVERVIEW 

One application where the use of two-dimensional autoregressive modeling has 
proven to be particularly advantageous is that of spectral estimation. The 2-D AR 
model was derived in Chapter II. This model is used to find the 2-D AR power 


spectral estimate P,(w;,w2) which is given by [Ref. 4] 
P,(w,wW2) = |H(u,w2)|? Pw , (5.1) 


where P,, is the power spectral density of the input and H(w ,,w2) is the transfer 
function of the 2-D AR model. The input is white noise with a constant power 
spectrum of amplitude o2. Therefore, we can write (5.1) as 


o2 


2 SS 
(wy Ww ) 1 als S (ke: ke €A) a(ky, kg )e~3 (wiki twek2)|2 


(5.2) 


The key to finding the 2-D AR power spectral estimate is determining the the pa- 
rameters a(k,,k2). The next section provides results of experimentation with the 
data-adaptive iterative algorithm as used to find the parameter estimates of a 2-D 
signal consisting of a sum of sinusoids in additive white noise. Results using the di- 
rect. inversion method and fixed data iterative Toeplitz approximation are provided 


for comparison. 


B. EXPERIMENTAL RESULTS 


The performance of the algorithms discussed in the previous chapters is com- 
pared by examining the estimated spectral densities computed from the parameters 


produced by each method. Each algorithm has been used to estimate the parameters 
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of the sinusoidal input array generated by 


x(n, 72) = Ay cos(27 fi1n1 + 2n fign2) 


+A. cos(27 foin1 + 2m fo2N2) + w(n, 5 n2) (5.3) 


where amplitudes A; = A2 = /2 and fi; are normalized frequencies in the range 
(0 < fi; < 0.5). The frequencies chosen are fi; = fi2 = 0.167 and fa = fo2 = 0.333. 
Additionally, the data-adaptive algorithm is tested with off-set frequencies to evaluate 
its peformance for that case. The noise term w(n1, 2) is zero mean and gaussian with 
the variance o? scaled give a desired signal-to-noise ratio (SNR). The SNR (in dB) 


is defined as [Ref. 4] 


2 


t 
2 ? 
Wasi 


N 
SNR = 10logiy >> 


7==] 


(5.4) 


where N is the number of sinusoids present and A is the amplitude term. The SNR 
is varied by holding A; constant and adjusting the value of o2. The algorithm’s 
performance using a 3 x 3 filter mask is evaluated for various sizes of input arrays 
and the common regions of support, including combined-quadrant (CQ). The results 
obtained using the data-adaptive algorithm are provided for the pre-windowed and 
covariance methods. The data-adaptive algorithm is tested with SNR’s of 10 dB and 
0 dB. It should be noted that the surface plots in all figures have been rotated 90 
degrees to show the separation of the spectral peaks more clearly. Contour plots are 
provided to better judge the accuracy of the estimates. The axes of the contour plots 
range from 0-30 which is taken from a 32 x 32 frequency domain array representing a 
quadrant of a 64 x 64 point 2-D FFT output. The 32 x 32 array covers the frequency 
range from 0 to z. 
1. Estimates by Direct Inversion 
The first case examined is that involving direct inversion of the autocorre- 


lation matrix to find the least squares AR parameters from the normal equations. An 


29 


8 x 8 input array was used with 1° quadrant quarter-plane support. ‘The resulting 


spectral density estimate is shown in Figure 5.1. This plot shows spectral peaks at 
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Figure 5.1: Spectral estimate using direct matrix inversion; SNR=10 dB. 


21.3 and 10.7 on the F1 and F2 axis respectively, which correspond to the normalized 
frquencies fi, = fig = 0.167 and f; = fez = 0.333 of the input data. The true 
locations of the frequencies are indicated by crosses. 
2. Spectral Estimates using the Iterative Method for Fixed Data 
a. Quarter-Plane Support | 

Figure 5.2 shows the spectral density of a(n,,n2) as estimated using 
the fixed dala iterative method with first quadrant QP support. As in the previous 
case Lhe frequencies estimated are very close lo the true frequencies. It is interesting 
to note however, that contrary to what might be expected, the quality of the estimate 
appears to degrade as the number of iterations is increased. The estimate produced 


after one iteration is clearly superior to the estimate obtained afler ten iterations. 
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This phenomenon is consistent with that experienced in the 1-D case whiere the 
iterative solution would begin to diverge after a certain number of iterations in cases 


with sinall dala sets [Ref. 12]. 
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Figure 5.2: Fixed Data, Ist Quadrant QP Support, SNR=10 dB. 


Results using second quadrant QP support proved to be generally 
better for this placement of the sinusoids than those obtained from the first quadrant. 
An example is provided in Figure 5.3. 

b. Combined Quadrant Support 

Spectral density estiinates using QP support demonstrate a tendency 

to distort elliptically. This distortion is evident in the first quadrant QP case shown 


in Figure 5.2. A method to compensate for this elongation of the spectral peaks 
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Figure 5.3: Fixed Data, 2nd Quadrant QP Support, SNR=10 dB. 
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involves combining the first and second quadrant estimates [Ref. 13] to give 


; 20° 
= ce eee oh 7 . (5.5) 
ge + P92 
The plots in Figure 5.4 indicate that this method yields accurate estimates. However, 


despite some improvement, the elliptical elongation at the base of the spectral peaks 


resulting from the influence of the first quadrant estimate is still quite noticeable. 
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Figure 5.4: Fixed Data, Combined Quadrant Support, SNR=10 dB. 


3. Spectral Estimates using the Data-Adaptive Iterative Method 
a. Pre-Windowed Input Arrays 
The input arrays were pre-windowed by adding P,; — 1 rows and 
columns of zeros to the top and left side of the array, where P; is the dimension 
of the filter mask. The purpose of pre-windowing the data is to compute the pa- 


rameter estimates using all the input data. Figure 5.5 shows the spectral estimates 
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for 1°¢, 274 and CQ support with an 8 x 8 input array. In Figure 5.6 we present the 





dad Quadrent 
fet Ovedreat 





Figure 5.5::8 x 8 input array (pre-windowed), SNR=10 dB. 


spectral estimates when only 16 data points (4 x 4 array) were available. The spectral 
estimates obtained from the 8 x 8 input array were quite acceptable, but some bias 
mevent in the first quadrant estimate which is carried over to the CQ estimate. 
The estimates obtained from the 4 x 4 array demonstrate the value of computing the 
CQ estimate. In this case the true signal frequencies were unrecognizable until the 
combined quadrant estimate was obtained. 

A final test of the pre-windowed algorithm consisted of running the 
filter mask only half way through an 8 x 8 input. The estimates are computed from 
the upper 4 x 8 section of the observed data. The result of this experiment is given in 


Figure 5.7. These results are similar to that of the 4 x 4 case and serve to reinforce the 
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Figure 5.6: 4 x 4 input array (pre-windowed), SNR=10 dB. 
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Figure 5.7: 4 x 8 portion of the input array (pre-windowed), SNR=10 dB. 
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inportance of the Combined Quadrant technique for 2-D AR spectrum estimation. 
b. Estimates using the Covariance Method 
As explained in section B of Chapter II the covariance method uses 
only the observed data. As a result, Py; — 1 rows and columns of the input array are 
not available when computing the parameter estimates. However, the experimental 
results obtained in this work indicate that the covariance method achieves better 
spectral peak resolution. Additionally, less bias is observed in the spectral estimates. 


This characteristic of the covariance method is evident in Fig 5.8. In this case the 
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Figure 5.8: 8 x 8 input array (Covariance Method), SNR=10 dB. 


best spectral estimate is clearly provided by the Combined Quadrant method. The 
bias present in the pre-windowed first quadrant support example has been completely 


eliminated by using the covariance method. 
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The next example demonstrates a unique characteristic of the data- 
adaptive iterative Toeplitz approximation algorithm. In this example the spectral 
estimate is computed using a 3 x 3 filter mask with a 4 x 4 input array. This means 
that only four elements of the input data are available to the filter to determine 
the parameter estimates. This result is of particular interest because the correlation 
matrix is actually rank deficient, and the problem cannot be solved by direct inversion. 
Nevertheless, the iterative Toeplitz approximation method does provide a means to 


compute the estimate. The results are given in Figure 5.9. Although the spectral 
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Figure 5.9: 4 x 4 input array (Covariance Method), SNR=10 dB. 


estimates for the first quadrant and Combined Quadrants are somewhat weak, the 
second quadrant estimate exhibits remarkable resolution and accuracy. 


As a further test of the algorithm, a signal with offset frequencies 


38 


fir = 0.250, fi2 = 0.125, fo: = 0.300, and f22 = 0.400 is generated and processed. 
The spectral estimates resulting from the computed model parameters are plotted 


in Figure 5.10. For this case, the best spectral estimate was obtained using first 
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Figure 5.10: 8 x 8 input array, offset frequencies, (Covariance Method), 
SNR=10 dB. 
quadrant suppport. 

The spectral estimates for a signal-to-noise ratio of 0 dB is ploted 
in Figure 5.11. The accuracy of the estimates is noticeably degraded in this case. 
However, a reasonable indication of the signal frequency content can be seen when 
the combined-quadrant method is used. 

The final case using NSIIP support is included for completeness. The 
filter mask had dimensions P, = P, = 4 and Lz = —2. The parameters were estimated 


for a 16 x 16 input array. The NSHP spectral estimate is plotted in Figure 5.12. The 
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Figure 5.11: 8 x 8 input array (Covariance Method), SNR=OdB. 





Figure 5.12: 16 x 16 input array (Covariance Method), SNR=10 dB; non- 
symmetric half-plane support. 
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resulting spectral estimate shown in Figure 5.12 is quite accurate, but this comes 
at the cost of greater complexity for the filter. NSHP support may be used as an 


alternative to combined-quadrant support in some cases. 


C. SUMMARY 


The data-adaptive iterative Toeplitz approximation algorithm has proven to 
be a viable means for 2-D AR spectrum estimation. Numerous factors must be 
considered when using this algorithm. In particular, the filter mask size should be 
chosen to correspond to the model order of the signal being processed. Although 
our overall experience has found that that CQ support produced the best results, 
no region of support seemed to have a clear advantage over the others for all cases. 
This indicates that the optimum region of support must be evaluated for the specific 


application that the algorithm is implemented for. 
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VI. RESTORATION OF IMAGES 


A. OVERVIEW 

The characteristics of typical images vary greatly from one region of the image 
to the next. For example, a region of an image containing a crowd or a building may 
have detailed variations in intensity, while another region representing the sky in the 
background will be essentially uniform in intensity. Two-dimensional data-adaptive 
filtering is particularly suited for processing data of this nature. The idea of adapting 
the processing to the local characteristics of an image is advantageous for many image 
processing applications, including image enhancement and restoration. 

There are two approaches to adaptive signal processing. The first approach 
involves adapting the filter output for each pixel. This is known as pizel-by-pirel 
processing [Ref. 4]. In this scheme, the processing method is based on the local 
characteristics of the image, degradation, and any other pertinent information con- 
tained in the pixel’s neighborhood region. This approach offers the greatest degree 
of flexiblity to the adaptive process, but has the highest computational complexity. 

When a more computationally efficient method is desired, a second approach 
known as subimage-by-subimage or block-by-block processing can be used [Ref. 4]. In 
this approach, the image array is divided into subimages and space-invariant tech- 
niques are used to process the data. This method can result in some discontinuities 
being present in the processed image. The size of the subimages directly affects the 
quality of the output image. The block-by-block method provides a faster means of 
processing, however, and may be desirable in some applications. 


In the next two sections of this chapter the data-adaptive Toeplitz approxima- 


42 


tion algorithm is implemented in image noise cancellation and line enhancer modes. 
For both of these applications the algorithm has been modified to perform pixel-by- 


pixel processing. 


B. IMAGE NOISE CANCELLATION 

The usual method of finding the estimate of a signal in noise is to pass the 
corrupted signal through a system that serves to suppress the noise while leaving 
the desired signal relatively unchanged. The noise canceler developed by Widrow 
(Ref. 14] is an example of such a system. Adaptive noise cancellation is a variation 
of optimal filtering that can be used for image restoration to some advantage. In 
Bee oalar when a recieved image is obscured by noise or another image transmitted 
from another location as shown in Figure 6.1, an adaptive noise canceler can be 


employed to extract the desired image from the composite signal. The noise canceler 


unwanted transmission 





desired transmission primary sensor 


Figure 6.1: A possible scenario for application of a noise canceler as used 
for image restoration. 


makes use of an auxiliary or reference input derived from one or more sensors located 
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at points in the noise field where the desired signal is weak or undetectable. This 
input is filtered and subtracted from the primary input containing both desired and 
unwanted signals. A block diagram of the noise canceler is shown in Figure 6.2. The 
reference signal and the undesired part of the primary input are correlated. For this 


reason, the unwanted signal or noise is eliminated or attenuated by cancellation. 


PRIMARY SENSOR 





INPUT 
» 
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INPUT 
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Figure 6.2: Adaptive Noise Canceler Block Diagram. 


To evaluate the performance of the data-adaptive iterative Toeplitz approxi- 
mation algorithm used in a noise canceler configuration, a 256 x 256 composite im- 
age array was created from two distinct images. This corrupted image is shown in 
Figure 6.3. The algorithm described in Chapter IV was modified to include a cross- 
correlation term given by | 


n—-1l 





Yr, = Ara-1 + ( ) daXn (6.1) 


where d,, is the element of the composite signal corresponding to the n** element of 
the reference signal being processed by the filter. This cross-correlation term is used 


to form the initial parameter estimate for each pixel being processed. 
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Figure 6.3: Corrupted image array provided as primary input to noise 
canceler. 


The resulting error e(n,,n2) computed at the output of the system becomes the 


restored image. This relationship is given by 
ARE n2) eae [d(nx, n2) 25 u(ny, n2)] ~~ u(ni, n2) ’ (6.2) 


where d is the desired image, u is the unwanted image or noise and w is the filtered 
estimate of uw. 

In the simulation the noise signal added to the desired image is scaled and 
lowpass filtered to simulate the effect of having the primary and reference sensors 
physically separated. The output of the noise canceler is shown in Figure 6.4. The 
unwanted image has been noticeably attenuated and the desired image is clearly 
visible. The original images are provided in Figure 6.5 as a reference for evaluation 


of the algorithm’s performance. 


C. ADAPTIVE LINE ENHANCER 


A special case of adaptive noise canceling occurs when only the signal that is 
corrupted by noise is available. In this case the recieved signal is delayed and used as 
the reference signal. The reference signal may be expressed as y(n,,n2) = r(n,;—6.n2), 


where 6 is the amount of spatial delay used. The block diagram for this system is 
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Figure 6.4: Processed image array after noise cancellation. 





Figure 6.5: The original images - undesired and desired. 
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shown in Figure 6.6. The main function of the delay parameter 6 is to remove any 





Figure 6.6: The adaptive line enhancer block diagram. 


correlation that may exist between the noise component in the original input signal 
and the noise component of the delayed adaptive filter input [Ref. 10]. The filter will 
respond by cancelling any components of the main signal z(n,,n2) that are in any 
way correlated with the secondary signal y(n;,n2) [Ref. 11]. The remaining noise 
component at the output of the filter is then subtracted from main signal resulting 
in the removal of the additive noise in the signal. In general, the delay 6 should 
be chosen such that it is approximately half the length of the correlation sequence 


(Ref. 11]. The input image and output of the adaptive line enhancer are shown in 


Figure 6.7. 


D. SUMMARY 
In this chapter, it has been demonstrated that the data adaptive Toeplitz ap- 
proximation algorithm may be used in image processing applications. While pixel-by- 


pixel processing was used in the above examples, it may be desirable to combine this 
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Figure 6.7: Input and output of the data-adaptive line enhancer. 


approach with block-by-block processing to obtain better results. This would entail 
dividing the image into smaller blocks before processing and then recombining the 
processed blocks. Additionally, reprocessing of the array estimates may be a viable 


means of improving image quality. 
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VII. CONCLUSIONS 


The data-adaptive iterative Toeplitz approximation algorithm is readily im- 
plemented from the iterative methods previously used for fixed data-arrays. The 
experimental results indicate that this algorithm is well-suited for a wide variety of 
2-D signal processing applications. Practically speaking, the iterative method can 
be successfully implemented for any application where the use of Wiener filter based 
adaptive filters has been successful. While more extensive evaluation is necessary 
before this method can be said to have a clear advantage over other methods, there 


is enough evidence of its capabilities to warrant further investigation in this area. 


A. PERFORMANCE EVALUATION SUMMARY 

As stated earlier, autoregressive model parameters may used to estimate the 
spectral content of 2-D arrays with high resolution. In experiments the performance 
of the data-adaptive iterative algorithm matched that of the direct inversion method, 
and in some cases, improved upon the performance of the fixed data iterative method. 
This improvement may be due to the way in which the iteration is carried out in the 
adaptive case versus the fixed case. In the fixed case the iteration is carried out using 
a correlation matrix formed from all the data at once. In the data-adaptive case the 
correlation matrix is formed recursively with new data incorporated at each itera- 
tion. Of particular note, was the ability of the data-adaptive algorithm to estimate 
frequencies using very small data sets and its ability to produce an estimate even 
when the correlation matrix is not of full rank . Additionally, its performance in a 
high noise environment was noteworthy. 


The algorithm also proved to be a viable method for data-adaptive image 
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restoration. The nature of this application differs greatly from that of spectral estima- 
tion in that the data sets are much larger and generally less correlated. Encouraging 
results were obtained for the data-adaptive noise canceler and line enhancer with only 


minor modifications of the same algorithm that was used for the spectrum estimates. 


B. RECOMMENDATIONS FOR FUTURE STUDY 


The primary objective of this thesis was to develop the data-adaptive iterative 
Toeplitz algorithm and obtain experimental results to evaluate its potential for use in 
applications that involve 2-D AR parameter estimation. There are many aspects re- 
garding this work that require more in-depth study and several applications that may 
benefit in the use of this algorithm. In particular, more detailed analysis of model 
orders and their relation to the algorithm’s ability to form spectral estimates from 
very small input arrays needs to be carried out. The use of singular value decompo- 
sition (SVD) and eigenvalue analysis of the observed data may provide some insights 
regarding the performance of the iterative algorithm. Also, further investigation as to 
why one region of support (ROS) provides greater resolution than another, may yield 
information that would permit the choice of an optimum ROS prior to processing 
the data. Further analysis for choosing values of the forgetting factor A and how it 
relates to the statistical nature of the data, is desirable as well. This may lead toa 
scheme for on-line adjustment of A, which would be particularly advantageous when 
processing images. Other future work could explore the use of different data win- 
dowing schemes and data-ordering schemes for passing the filter mask over the data. 
One possible method that has not been investigated is an expanding square starting 
at one corner of the array. This method seems well suited for block-by-block image 
processing. Finally, the data-adaptive iterative Toeplitz approximation algorithm has 


the potential for use in non-linear applications, such as the identification of Volterra 
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